arXiv:1504.02399v2 [cond-mat.quant-gas] 20 Aug 2015 


Time-of-flight images of Mott insulators in the Hofstadter-Bose-Hubbard model 

M. Iskin 

Department of Physics, Kog University, Rumelifeneri Yolu, 34450 Sariyer, Istanbul, Turkey. 

(Dated: August 21, 2015) 

We analyze the momentum distribution function and its artificial-gauge-field dependence for the Mott insu¬ 
lator phases of the Hofstadter-Bose-Hubbard model. By benchmarking the results of the random-phase approx¬ 
imation (RPA) approach against those of the strong-coupling expansion (SCE) for the Landau and symmetric 
gauges, we find pronounced corrections to the former results, which is a clear manifestation of the critical role 
played by quantum fluctuations in two dimensions. 
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Introduction: The momentum distribution function n(k) 
of atoms, which is defined as the Fourier transform of the 
one-body density matrix, can be directly measured in cold- 
atom systems by time-of-flight absorption imaging of freely 
expanding gas (ll-l^. Since these systems are extremely di¬ 
lute, the atom-atom interactions are negligible during such an 
expansion, and the position of atoms at time r are strongly 
correlated with their velocity distribution at the moment of 
release from the trap, i.e., r = fikr/m with H the Planck con¬ 
stant and m the atomic mass. Therefore, the n(k) of atoms 
has not only been the easiest observable to measure but also 
been routinely used for probing distinct phases of matter in 
atomic systems. 

In addition, followed by the recent advances in creating ar¬ 
tificial gauge fields in atomic systems flUt], there has been 
growing interest in first the realization of the Hofstadter-type 
lattice Hamiltonians and then the detection of the resultant 
many-body phases IMl- For instance, the MIT group has 
in their latest preprint measured the n(k) of atoms in the su¬ 
perfluid (SF) phase im, revealing both the reduced symme¬ 
try of their specific gau ge held and the resultant degeneracy 
of the ground state iflin . There is no doubt that such a ca¬ 
pacity to tune strong gauge fields together with strong inter¬ 
actions paves ultimately the way for creating and observing 
uncharted many-body phases and transitions in between, one 
of the immediate candidates of which is the renowned SF-MI 
transition iQli. 

Motivated by these recent works, in this brief paper, we 
study n(k) of atoms for the MI phases of the Hofstadter-Bose- 
Hubbard model on a square lattice. For this purpose, we com¬ 
pare the results of RPA and SCE approaches for the Landau 
and symmetric gauges, and And substantial corrections to the 
former results depending strongly on the specified gauge. 

Hamiltonian and Phase Diagram: These results are ob¬ 
tained for the following Hamiltonian 

H = — ^ ^ ^ ^ ni(ni — 1) — p 'y ^ m, (1) 

ij i i 

where the hopping parameter Uj = connects nearest- 

neighbor sites with phase factor 0^ taking the gauge fields 
into account, c| (ci) creates (annihilates) a boson on site i, 
the boson-boson interaction is on-site and repulsive U > 0, 
rii = cl Ci is the number operator, and /r > 0 is the chemical 
potential. In this paper, we compare the results of the usual 


(A) no-gauge limit, where Oij = 0 for all hoppings; with those 
of {B) Landau gauge, where Oij = 2Tt(l)u for (u, v) to (u, v + 
1) and 0 for {u,v) to (u + l,u) hoppings; (C) symmetric 
gauge, where Oij = irtpu for {u,v) to {u,v + 1) and —ndv 
for {u,v) to {u + l,u) hoppings; and (D) MIT gauge lllffl . 
where Oij = 2TT(j>{u + v) for {u,v) to {u,v + 1) and 0 for 
(u, t;) to (u + 1, v) hoppings. Here, (u, v) corresponds to the 
Cartesian coordinates of site i, and 0^ are chosen such that 
the magnetic flux tf) = p/q is the same for all gauges, where p 
and q are co-prime numbers with p < q. 



FIG. 1. (Color online) The ground-state phase diagram is shown as a 
function of chemical potential p,, magnetic flux ()) = p/q and hopping 
strength 4f. 


In the atomic {t = 0) limit, since H commutes with fii, 
the thermal average rii = (ni) is such that the ground-state 
energy is minimised for a given p, leading to a uniform oc¬ 
cupation (rii = n) of bosons thanks to the translational in¬ 
variance of H. When U — 0 and p = 0, the spectrum of H 
corresponds to the celebrated Hofstadter butterfly |fl^[T^ . It 
is also very well-known that the range of p about which the 
ground state is a MI with an integer occupation n decreases 
as a function of increasing t/U, and depending on n and (j), 
the Mis disappear at a critical value of t/U, beyond which 
the system becomes a SF Q. For instance, the qualitative 
phase diagram of H can be obtained within the mean-field 
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approximation, e.g., the decoupling or variational Gutzwiller 
techniques, leading to |fl5l - [l^ 


1 n + 1 n 

ePt Un — U{n — 1) — fj. 


( 2 ) 


at zero temperature for the MI-SF phase transition boundary, 
where n > 0 is an integer number. Here, is the mini¬ 
mal eigenvalue of the hopping matrix 
and it corresponds to the maximal single-particle kinetic en¬ 
ergy of the Hofstadter butterfly, e.g., e° = At when ^ = 0. 
Since the effects of 9ij enter Eq. (|2| through its dependence 
on eP'^, the mean-field phase boundary is clearly independent 
of the gauge, which is simply because only the position in the 
magnetic Brillouin zone but not the value of eP'^ depends on 
the gauge. However, this is not the case for the SF properties 
which are gauge dependent within the mean-field approaches. 

In Fig. [T] we show the ground-state phase diagram as a 
function of n, <p = p/q and At, which is obtained by solving 
Eq. (|2l) together with the Harper’s equation. Both the sym¬ 
metry around p/q = 1/2 and the intriguing structure of the 
MI-SF phase transition boundary are due to the dependence 
of eP9 on IHIIl. In addition, the incompressible (com¬ 
pressible) MI (SF) phase grows (shrinks) when (f) increases 
from 0, a consequence of which is due to the localizing ef¬ 
fects of magnetic flux on particIeSjmd all of these results are 
in agreement with earlier findings ■ Having introduced 

the model Hamiltonian and reviewed its phase diagram, next 
we are ready to discuss the momentum distribution of bosons 
for the Mis. 

Momentum Distribution: As discussed in the Introduction, 
the n(k) of atoms corresponds to the Fourier transform of the 
one-body density matrix, and it is given by ifTsl - l^ 


i(k) = 


M 


33 ' 




(3) 


where M is the number of sites and rj = {ua, va) is the po¬ 
sition of site j with a the lattice spacing. In the following, we 
set the Fourier transform of the Wannier function r(;(k) to 1, 
since it depends on the particular optical lattice potential and 
has nothing to do with our H. 

In this paper, we calculate n(k) for the Mis using two ap¬ 
proaches: (/) RPA JUm and {II) SCE in t/U ifllfll. 
We emphasize that while the result of the RPA approach cor¬ 
responds to the exact n(k) only in the limit of infinite dimen¬ 
sions and zero magnetic flux, the results of the SCE approach 
are exact in two dimensions for the specified gauges up to the 
given order in t/U. 

(I) Random-Phase Approximation: In the RPA ap¬ 
proach in Si, since the thermal averages of products of 
operators are replaced by the product of their thermal aver¬ 
ages, the fluctuations are not fully taken into account. After a 


lengthy but straightforward algebra, one finds 


'■RPA 


1 


9-1 




£r(k) + U 

^[ef(k)]" + 2(7£f(k) + C/2 



for a MI with n bosons per site at zero temperature, where 
U = U{2n A- 1) and £f'^(k) is the energy dispersion of a sin¬ 
gle particle in the fth Hofstadter band. Note that the form of 
Eq. dUi is exactly the same as the usual Bose-Hubbard model, 
i.e., the main difference is a sum over the Hofstadter bands, 
and that it has an overall factor of 1 /q in comparison to the 
one given in Ref. 1^ . While the set of £^'^(k) values de¬ 
pends only on (/> and lattice geometry, their corresponding po¬ 
sitions in the 1st magnetic Brillouin zone, and therefore n(k), 
are gauge dependent H^IU. For instance, n(k) exhibits q 
peaks as a function of k, and only the number q but not the 
positions are controlled by (j). Note that eP‘^ = max{£^‘^(k)} 
in Eq. (I2]i which is also a gauge-independent quantity as re¬ 
marked above. In particular, when (/ = 0, a d-dimensional 
hypercubic lattice gives rise to a single band with dispersion 
£°(k) = —2f cos(fcia), and it is already established 

that r ip PA (k) becomes exact as d —>■ oo while keeping dt 
fixed illdl- 

To compare Eq. (|4| with our exact results of the SCE ap¬ 
proach derived below, let us expand n^p^(k) in a power series 
up to 3rd order in t/U, leading to 


^ e=o 

+ 3n(n+l)(2n+l) g^^p,(k)]2 
^ e=o 

An{n -\- l)(5n^ -f 5n -f 1) ^ ,^.|3 

^ e—f\ 


(5) 


For a given </>, the sums over Hofstadter bands can be eas¬ 
ily evaluated for a given gauge by noting = 

Trace{ [T^’‘?(k)]* }, where T^’‘?(k) describes the kinetic en¬ 
ergy of a single particle in the 1st magnetic Brillouin zone. 

For instance, TP'J(k) is a q x g matrix in the Landau 
gauge EQ 

do Q ^ _^g-zfex,a 

di . 0 

— tei-k^a Q _ _^g-ifcxO 

(6) 

with hg = —2tcos{kya -f 2TT(j)i), for which the s = 1 
trace equals to —2f[cos(fca;a) -f cos(fcya)] when {p, q) = 
(1,1) and it vanishes for q > 1; the s = 2 trace 

equals to [cos(fca,a) -I- cos(fcya)] when (p, q) = (1,1), 
to [cos^(fca;a)-f cos^(fcya)] when (p, q) = (1,2) and 
to 4qf^ for q > 2; and lastly the s = 3 trace equals to 
—[cos(fca;a) + cos(fcpa)]^ when (p, q) = (1,1) and to 
—6t^ [cos(3fcxa) + cos(3fcya)] when q = 3, but it vanishes 
for q > 3. Thus, Eq. (|5]l reduces to 












3 


'^RPA 

(k) 

= n 




32n(' 

^^RPA 

(k) 

= n 

^ ( 

p3 

^RPA 

(k) 

= n 

-h - 

,P,g>3 

RPA 

(k) 

= n 

-h - 


4n(n + 1) 

iu/t) 


cos{kxa) + cos(fcya)] + 


12n{n + l)( 2 n + 1 ) 


{U/tf 

6n{n + l)( 2 n + 1 ) 


(u/ty 

cos{kxa) + cos(/cya)]^, 


cos{kxa) + cos{kya)Y 


{U/tf 


cos{2kxa) + cos{2kya) + 2] + 0{t/Uy, 

cos{3kxa) + cos(3fcyo)], 


(U/ty 

12n(n + l)(2n + 1) 8n{n + l)(5n^ + Sn + 1) 


{u/ty 


(u/ty 


(7) 

( 8 ) 
(9) 

( 10 ) 


Equations (fTlfTOl) clearly show that the first k dependence of 
n^p^(k) arises at the gth order in t/U. More importantly, we 
note that Eqs. (ITIfTOl) are symmetric in kx and ky even though 
the spatial symmetry between x and y directions is explicitly 
broken by the Landau gauge. Note also that Eq. (O coincides 
with that of the cj) = 0 result since (k) is a periodic func¬ 
tion of (j) with a period of 1 dB. Unlike the tp = 0 case 
for which the RPA approach captures the essential features of 
n°(k) even in finite dimensions dlii, next we use the SCE 
approach and show that the corrections to npp^(k) are quite 
dramatic in the presence of gauge fields in two dimensions. 

(II) Strong-Coupling Expansion: In the SCE approach id 
B, the wave function of Mis is achieved via a many-body 
perturbation theory in the kinetic energy term up to 3rd order 
int/U. In principle, one can apply the perturbation theory on 

the Oth-order wave function I'I'mI) = O^i |0)/v^, 

where |0) is the vacuum state, and calculate I'I'mi) up to any 
desired order. However, since the number of intermediate 
states increases dramatically, here we perform this expansion 
only up to 3rd order in t/U, and obtain I^/mi) = \'4’mi)/A 
where 




m'm" m" 


m mi 

\^rn"'m" ^m"m' ^m'O 
^Om'" EoYn,/ 


/\ , ^m"m'^m'0 \ //\ 

-1^ ) 


Eom"Eom' 


K") + ' 


( 11 ) 


is the unnormalized wave function which needs to be di¬ 
vided by a proper normalization coefficient A in order 
to get the correct order of perturbation. Here, Vm'o — 
— yYjjj' 1 ' 1 'mi) connects the Ist-order intermedi¬ 
ate states |m') to is their Oth- 

order energy difference, and \m'') and \m"') are respectively 
the 2nd and 3rd-order intermediate states. Note that while 
and \m'), |to') and |to"), and \m") and \m"') states 
are connected to each other with a single hopping, \m") and 
\m!") states must be different from the state. There¬ 

fore, the normalization condition (4 'miI^'mi) = 1 gives = 
1 -I- 4n(n -I- l)Mt^/U^ 0(t/U)‘^, which has vanishing 1st 

and 3rd order terms. 


After a very lengthy and tedious algebra, one finds 


/ T it I -r \ r 2n{n 1) 
(WMi|a],aj|VMi) = ndjf H- — - tjj> 


-h 


3n(n -I- l)(2n + 1) 
IP 


^331 ^3. 


31 


4n(n -I- l)(5n^ -I- 5n -|- 1) 


^332 ^3231 ^3 


3132 


n{n -I- l)(131n^ -I- 131n -|- 26) 
IP 


t tjji 


( 12 ) 


for a square lattice with nearest-neighbor hopping at zero 
temperature. We note in Eq. (fT^ that the 2 terms that are 
explicitly proportional to P are finite-d corrections, includ¬ 
ing the 2nd term in the 2nd line and the 4th line, as they 
vanish in the d —t oo limit while keeping dt fixed. Since 
Eq. (IT^ is derived exactly using a generic hopping matrix 
tij, we are ready to benchmark it against the results of the 
RPA approach for a number of specified gauges. Eor this 
purpose, we make use of the following identities: the sum 
Y^eZo cos(a — 2mr(j)£) equals to qcos{a) when q = n and it 
vanishes for q > w, the sum YYo cos^(a — 2'K(/>t) equals to 
gcos^(a) when (p, g) = {(1,1), (1, 2)} and g/2 for g > 2; 
and the sum YYo cos^(a — 2'K(j)t) equals to cos^(q;) when 
(p, g) = (1,1) and to 3 cos(3q;)/4 when q = 3, but it vanishes 
for(p,g) = (1,2) org > 3. 

(II-A) No-Gauge Limit: Setting = 0 for all hoppings in 
Eq. (O, we obtain 


On t 2n(n-h 1) on \ 

n (k) = n- — - 1 (k) 

3n(?T.-|-l)(2n-|-1) ^^2 


C/2 

4n(n -|- l)(5n2 q- 571 -|- 1 ) 

JP 

n(n-I-l)(131n^-I-131n-I-26)^2 0 

IP 


{[eyk)y- 4 y} 

e°(k)]3 

(13) 


where e°(k) = —2t [cos(fca;a) -I- cos(A:ya)] is the usual dis¬ 
persion relation for a square lattice. Since the two terms that 
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are explicitly proportional to are finite-d corrections, they 
are not captured by the result of the RPA approach that is given 
in Eq. (|2l). 


(II-B) Landau Gauge: On the other hand, setting 9ij = 
2'iT(j)u for {u,v) to (u,u + 1) and 0 for (m,u) to [u + l,w) 
hoppings in Eq. (fT^ . we obtain 


11 4n(n+l) 12n(n + l)(2n +1) 


Ul (k) = n + 

+ 


{U/t) 

32n(n + l)(5n^ + 5n + 1) 

WNf 


cos{kxa) + cos{kya)] + 


(u/ty 


{[cos(fca;a) + cos{kya)y — l} 


(14) 


[cOs(A:a:a) + cos(fcya)]'^ — 


3 2n(n + l)(131n^ + 131n + 26) 


i 2 i, 1 ^ 4n(n + 1) ^ 6 n(n + l)(2n + 1) 


Ul (k) = n + 


{U/t) 


■ cos{kxa) 


{U/ty 


{U/t)^ 
cos{2kxa) + cos{2kya)] 


cos{kxa) + cos(fcya)], 


(15) 


+ 


32n(n + l)(5n^ + 5n + 1) 


{U/ty 

^ ^ + {U/t) 


cos(fc2:a)[cos^(A;a;a) + cos^(fcya)] — 


2n{n + l)(131n^ + 131n + 26) 

ww 


cos{kxa), 


cos{kxa) + cos{2kxa) 

\p/u 


(16) 


8n(n + l)(5n^ + 5n + 1) /m \ r- n \^ 2n(n + l)(131n^ + 131n + 26) 

H--[cos(3Ka;a) + cos{3kya) + 6 cos(K 2 :a)J- nr u\?. - cos{kxa), 


{U/ty 


{U/ty 


P,q> 3 f.^ _, 4u(n + 1 ) ^ 6 n(n + l)( 2 n + 1 ) 


(k) = n + 


■ cos{kxa) 


+ 


{U/t) 

8n{n + l)(5n^ + 5n + 1) 

W/t) 


cos{2kxa) 


(17) 


{U/ty 

r /o 7 ^ r^T r, /n / M n \i 2n(n + l)(131n^ + 131n + 26) 

3 |cos(3fca;a) + [7 + 2cos(27rp/g)J cos(fc:ca)|- {U/t)^ - cos{kxa). 


Note that Eq. (fT4l i exactly coincides with Eq. (fOT l since cj) = 1 
and 0 are equivalent in this gauge. We also note that, unlike 
the results of the RPA approach that are given in Eqs. dlloli, 
these exact results are not symmetric in kx and ky, showing 
that it is only the first ky dependence that arises at the qth 
order in t/U. This is not surprising because while the one- 
body correlation operator cj,Cj connects IT'mi) to itself at the 
1 st order in x direction, the connection is established at the 
gth order in y direction due to the presence of 2'K<j)u. In ad¬ 
dition, on top of the RPA contributions, Eqs. (I14H17I ) contain 
various other terms, showing that the hnite-d corrections are 
quite substantial in the presence of gauge helds in two dimen¬ 
sions El. Thus, one of our main conclusions in this paper 
is that the mismatch between the results of RPA and SCE ap¬ 
proaches grows so dramatically as q increases from 1 that the 
former approach fails to reproduce any of the exact terms up 
to 3rd order in t/U for q > 3. 

(II-C) Symmetric Gauge: Similarly, setting 0ij = tk/u for 
(u, v) to (u, u-f 1 ) and —TT(j)v for {u, v) to (u-f 1 , v) hoppings 
in Eq. (fTSl l. we obtain 

11 12n(n-I-l)(2n-I-1) 1 

ns (k) = n -f--^{[cos(fc^a) 

+ cos{kya)y - 1} + 0{t/U)‘’', (18) 

nf’«>i(k) =n-f C>(t/[/)^. (19) 

Note that Eq. ( fTSb does not reproduce Eq. ( fT3] ) since </> = 1 


and 0 are not equivalent in this gauge. We also note that, un¬ 
like the results of the SCE approach for the Landau gauge that 
are given in Eqs. (I14H171 I. here the k dependence is not only 
symmetric in kx and ky, thanks to the spatial symmetry be¬ 
tween X and y directions, but also the first k dependence arises 
at the 2q\h order int/U. This is also not surprising because the 
one-body correlation operator c|,Cj connects IT'mi) to itself at 
the 2qth order in both directions due to the presence of ttc/u. 
In addition, the k-independent 2nd order term in Eq. (18) is 
a hnite-d correction to the result of the RPA approach in this 
gauge. Therefore, n 5 ^(k) becomes more and more featureless 
function of k as g increases from 1, especially deep in the Mis 
when t/U is very small. 

(II-D) MIT Gauge: Lastly, setting 9ij = 27r^(u + v) for 
{u,v) to {u,v + 1 ) and 0 for {u,v) to {u + l,u) hoppings 
in Eg . (fl^ leads exactly to Eqs. (I14II17I) . and therefore, the 
MIT iJol and Landau gauges have exactly the same n(k). 

Conclusions: To summarize, we studied the expansion im¬ 
ages of atoms for the MI phases of the Hofstadter-Bose- 
Hubbard model on a square lattice. In particular, we explicitly 
calculated the momentum distribution function for the Landau 
and symmetric gauges with both RPA and SCE approaches, 
and found marked corrections to the former results depending 
strongly on the specihed gauge. Such a comparison clearly 
manifests the importance of the critical role played by quan¬ 
tum fluctuations in two dimensions. 
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